function detectPower = evalSeq(transSeq,SOA, HRF_file, HRF_res)
%evalSeq: Calculate the detection power  of a given sequence
%
%
%   Inputs
%   -------
%   transSeq: Sequence of transition distances (-1 corresponds do NA
%       transitions)
%   SOA: Stimulus onset asynchrony, or inter-trial interval (in ms)
%   HRF_file: Path to the hrf .txt file
%   HRF_res: Temporal resolution of the hrf file (in ms)
%   
%   PS: NA distances in the transition sequence should be represented with -1
%
%   _______________________________________________
%   by Marcelo G Mattar (13/07/2012)
%   mattar@sas.upenn.edu

if nargin<4
    HRF_res = 100;
end

if nargin<3
    HRF_file = 'Not Entered';
end

if strcmp(HRF_file, 'Not Entered')
    rawHRF = [5.55112e-17,0.00767224,0.0162726,0.025849,0.0364457,0.0481024,0.0608536,0.0747281,0.0897482,0.10593,0.123281,0.141801,0.161484,0.182312,0.204261,0.227297,0.251376,0.276448,0.30245,0.329315,0.356963,0.38531,0.414261,0.443716,0.473569 ,0.503705,0.534007,0.564352,0.594614,0.624665,0.654374,0.68361,0.712242,0.740139,0.767176,0.793226,0.818169,0.841889,0.864277,0.885229,0.90465,0.922451,0.938554,0.952889,0.965396,0.976025,0.984737,0.991502,0.996304,0.999136,1,0.998913,0.995897,0.990989,0.984232,0.975681,0.965398,0.953452,0.939923,0.924892,0.908452 ,0.890696,0.871725,0.851641,0.830549,0.808556,0.78577,0.7623,0.738253,0.713734,0.688848,0.663695,0.638375,0.61298,0.5876,0.562319,0.537218,0.51237,0.487844,0.463702,0.44,0.41679,0.394115,0.372016,0.350523,0.329666,0.309467,0.289942,0.271104,0.25296,0.235516,0.218771,0.202721,0.187361,0.172683,0.158674,0.145324,0.132616,0.120536,0.109067,0.0981926,0.0878945,0.0781552,0.0689569,0.0602817,0.0521123,0.0444312,0.0372216,0.0304668,0.0241507,0.0182571,0.0127706,0.00767572,0.00295715,-0.00140033,-0.00541187,-0.00909271,-0.0124583,-0.0155242,-0.0183063,-0.0208208,-0.0230841,-0.025113,-0.0269245,-0.0285358,-0.0299644,-0.0312277,-0.0323428,-0.0333269,-0.0341965,-0.0349674,-0.0356548,-0.0362724,-0.0368327,-0.0373465,-0.0378227,-0.0382678,-0.038686,-0.0390787,-0.0394442,-0.0397777,-0.0400708,-0.0403115,-0.0404842,-0.0405693,-0.0405433,-0.0403786,-0.0400439,-0.039504,-0.0387199,-0.0376493,-0.0362464,-0.0344627,-0.0322468,-0.0295455,-0.0263035,-0.0224645,-0.0179715,-0.0127673,-0.00679528]';
else
    rawHRF = dlmread(HRF_file);
end


%% Mean center
% Mean center non-NA values
transSeq_meancentered = transSeq - ((transSeq~=-1) * sum(transSeq(transSeq~=-1))/sum(transSeq~=-1));

% Set NA values to zero
transSeq_meancentered(transSeq==-1) = 0;


%% Upsample transSeq
factor = SOA/HRF_res;
transSeq_upsampled = zeros(length(transSeq_meancentered)*factor,1);
for i=1:length(transSeq_upsampled)
    if mod(i-1,factor)==0
        transSeq_upsampled(i) = transSeq_meancentered((i-1)/factor + 1);
    else
        if ceil((i-1)/factor + 1) < length(transSeq_meancentered)
            transSeq_upsampled(i) = transSeq_meancentered(floor((i-1)/factor + 1)) + ((transSeq_meancentered(ceil((i-1)/factor + 1)) - transSeq_meancentered(floor((i-1)/factor + 1)))/factor) * (mod(i-1,factor));
        end
    end
end

%% Convolve transSeq with the HRF

% FFT both hrf and transSeq
hrf_fft = fft(rawHRF,length(transSeq_upsampled));
transSeq_fft = fft(transSeq_upsampled);

% Normalize hrf_fft
hrf_fft = hrf_fft/max(abs(hrf_fft));

% Apply notch filter
hrf_fft(f_axis(length(hrf_fft),0.1)<0.01) = 0;

% Multiply (convolve) transSeq_fft and hrf_fft
convolved_fft = hrf_fft .* transSeq_fft;

% IFFT the product
convolved = ifft(convolved_fft);

%% Calculate Detection Power
detectPower = var(convolved)/var(transSeq_upsampled);

%==========================================================================





function freq_axis = f_axis(numPoints, sampleRate)
    N = numPoints; 

    % Sampling frequency
    Fs = 1/sampleRate;

    % Calculate the number of unique points
    NumUniquePts = ceil((N+1)/2);

    % Create the first half of the frequency axis
    f_1sthalf = (0:NumUniquePts-1)*Fs/N;

    % Create the second half of the frequency axis (reverse order)
    if rem(N,2)
        f_2ndhalf = f_1sthalf(2:end);
    else
        f_2ndhalf = f_1sthalf(2:(end-1));
    end

    % Attach both halfs together to form the full frequency axis
    freq_axis = [f_1sthalf fliplr(f_2ndhalf)];